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We compute analytically the anisotropic flow in an expanding mixture of several species of rel- 
ativistic massive particles. We find that a single collision per particle in average already leads to 
sizable elliptic flow, with mass ordering between the species. 
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I. INTRODUCTION 

The hundreds or even thousands of particles in the 
final state of a nucleus-nucleus collision at high en- 
ergy are emitted asymmetrically in the plane transverse 
to the collision axis. This phenomenon, referred to 
as anisotropic flow, is conveniently quantified by intro- 
ducing the Fourier expansion of the (invariant) particle 
yield [l| . Restricting ourselves to the transverse momen- 
tum distribution, 
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where the Fourier coefficients are given by 
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In Eqs. ([T]) ip denotes the azimuth of outgoing particles 
with respect to the direction of the impact parameter. 
Anisotropic flow is thus a collective effect, that signals a 
correlation of each particle to the reaction plane. 

Given this collectivity property, a natural modeling of 
anisotropic flow consists in considering the particles as a 
continuous medium that evolves according to the laws of 
relativistic fluid dynamics 0] . This description implicitly 
assumes that particles undergo many collisions, to ensure 
(quasi-)local equilibrium. 

On the opposite, one may also investigate anisotropic 
flow far from equilibrium, when outgoing particles inter- 
act seldom, which is our purpose in this article. A small 
number of interactions per particle is expected in several 
cases. The most obvious one is that of peripheral col- 
lisions, where only relatively few particles are emitted: 
then the number of their reinteractions, integrated over 
the whole expansion, will be small. Other instances of 
presumably far-from-equilibrium evolution are the early 
and late stages of the expansion in collisions with arbi- 
trary centrality: either before the systems equilibrates, 
as was considered in Ref. 0, Q (here, the rate of col- 
lisions remains small for a short while, before becoming 
large); or at the very end of the evolution, when collisions 
become seldom around the kinetic freeze-out. 



In this article, we describe the expanding system as a 
mixture of several relativistic massive components gov- 
erned by the Boltzmann equation, which we solve analyt- 
ically in the limit of a small number of collisions per par- 
ticle, focusing on calculating the flow harmonics (jlbl) . To 
our knowledge, the only previous attempt in that direc- 
tion is that of Ref. @ , which however only deals with non- 
relativistic particles in the low-eccentricity limit. Here 
our results will hold across all centralities, and we are 
not restricted to low momenta. 

In Section [TT] we introduce our model and outline the 
computation of the average number of collisions per par- 
ticle and of the resulting anisotropic flow. Section IIIII 
contains the technical details of our calculations in the 
case of massive particles — the limiting cases involving 
massless particles are dealt with in Appendices [Bl and [Cl 
We present our results in Section |IV] and discuss their 
meaning, as well as the limitations of our model in Sec- 
tion [V] Eventually, Appendix lAl introduces for the sake 
of reference some of the properties of the special function 
in terms of which we express our analytical results. 



II. MODEL 

Our purpose in the following is to describe the evo- 
lution of the expanding system created in a high-energy 
nuclear collision in the limit of a small average number of 
collisions per outgoing particle. For that, we model the 
system as a two-dimensional mixture of various compo- 
nents with respective particle masses m;, which can in- 
teract elastically with each other. As we are interested in 
the anisotropics in the transverse expansion, we restrict 
ourselves to the transverse (x, y)-plane, neglecting the 
influence of the longitudinal expansion. This is a valid 
ansatz as long as we consider the qualitative aspects of 
our final results only. Our system will be "homogeneous" , 
i.e. when we consider a thermal-like momentum distribu- 
tion with inverse slope parameter T, then the latter will 
be the same everywhere. Eventually, we assume that 
the system is invariant under the parity transformation 
x — > —x in position space, and that this symmetry holds 
during the whole evolution. 

A suitable approach consists in using a Boltzmann de- 
scription, approximating its exact solution as a small 
modification to the solution of the collisionless equation. 
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To each component of the mixture we associate its dis- 
tribution function fi(t, x, p,), normalized to the particle 
number A/j. The corresponding particle density is 



m (i,x) 



- Jd 2 pifi(t,x, p^, 



and the (transverse) momentum distribution 
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Multiplying this momentum distribution by cos nipi, with 
ifi the azimuth of Pi, and averaging over ipi, one obtains 
the flow coefficient v n (pi,t) at time t [Eq. f|l"b[)] . The 
"usual" v n (j>i) is the large-time limit of v n (j>i,t). 

The distribution function /j = /j(t, x, p.^) for each com- 
ponent obeys the classical relativistic Boltzmann equa- 
fiord 
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With /fc = fk(t,X,p k ), f! EE /i(t,X,p9, f' k EE /fc(t,X,p' fc ), 

while the factor 1 — ^(5^ ensures that the equation holds 
for both identical and nonidentical particles. The relative 
velocity is given by 
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and der^/d© is a differential scattering cross section for 
collisions between particles of species i and k. Hereafter 
we shall consider only elastic collisions, so that the num- 
ber of particles for each species remains constant, with 
an isotropic, constant differential cross section, which we 
shall for the sake of brevity denote <7d- Note that in 
our two-dimensional approach the cross section has the 
dimension of a length. 

Integrating the Boltzmann equation over x yields for 
the first term in the left-hand side the time derivative 
of the momentum distribution @. On the other hand, 
given the assumed parity of the system, the spatial gra- 
dient V x /i is odd under the parity transformation, and 
its integral over x vanishes. The integrated Boltzmann 
equation thus describes the time evolution of the momen- 
tum distribution, in particular of its anisotropics, under 
the influence of elastic collisions 

In the absence of rescatterings, the Boltzmann equa- 
tion is satisfied by the free-streaming solutions 



/f ) (t,x,p J ) = /f ) (0,x~v l t,p l ), 



(6) 



which only depend on the initial distribution at t = 0. 
If the latter is isotropic in momentum space, the free- 
streaming solution remains isotropic, i.e. it does not de- 
velop any anisotropic flow. As initial condition, we shall 
consider a factorized distribution, with a Gaussian de- 
pendence in position space and an isotropic distribution 
fi(pi) in momentum space: 
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with pi = | pi | and the normalization JdpiPifi(pi) = 1, 
and R y > R x , i.e. the reaction plane lies along the x-axis. 
Note that this distribution is invariant under the x — > — x 
transformation, and that with the assumed non-parity- 
violating cross section the property holds throughout the 
evolution. 

Hereafter we shall rewrite the mean square radii as 
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that is 
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e is thus the usual eccentricity of the initial distribution. 

When collisions are present, the distribution function, 
starting with the same initial condition at t = 0, changes 
to 



fS^Pi) = ri 0) (t,x,pi) + /,Hi,x,p 2 ) + 



with /■ ^ much smaller than /v 



(9) 



Due to the colli- 
sions, fi may develop anisotropies in momentum space, 
which, neglecting the next terms in the expansion, will 
be those of / 4 . Reinserting expansion (|9]) in the Boltz- 
mann equation, one should as a first approximation con- 
sider only fj ' + in the left-hand side and in 
the right-hand side. Our approach thus differs from the 
traditional Chapman-Enskog formulation, in which the 
collision integral of the leading term in the expansion 
vanishes 0, Chapter 6]. 

Let us now focus on the momentum anisotropies. The 
gain term in the collision integral does not contribute 
to developing such anisotropies, at least to the assumed 
order of approximation, thanks to the isotropy of the 
cross section. More precisely, the gain term involves the 
free- streaming particle distributions after the collision. 
The "memory" of the latter only extends back to the 
time of their last collision — technically, this amounts 
to replacing t = (resp. t) by the collision time i co n 
(resp. by t — t co \\) in the right-hand side of Eq. © - 



We neglect any long-range interaction between particles, which 
could be implemented as a mean-field term, as considered in the 
context of anisotropic flow in Ref. [y| . 



2 The "small parameter" that controls the expansion is the average 
number of collisions per particle, which is proportional to a^. 



3 



and thus bears no correlation to the initial geometry — 
which means that the product f[ does not involve 

the azimuths ifi and <fk of the momenta, so that the 
integral of fi f^' v ik cosnipi over both azimuths gives 
0. 

To derive the numerator of Eq. (jlb|) . we thus just have 
to consider the loss term, proportional to —fj f^Vikcr^. 
Let us first not perform the integration over p&: the term 
then depends only on the velocities v^, of the collid- 
ing particles. Multiplying it by cos nipi, we can then in- 
tegrate over time and the position, as well as over the 
azimuths (fk and tfi . This yields a function, hereafter re- 
ferred to as the "anisotropy of the velocity distribution" , 
which integrated over |pfe| yields the numerator in the 
definition (llb|) of the flow harmonic v n {pi) for particle i. 

Anticipating on one of our results, this distribution 
involves the product of fi(pi)f k {pk) with a function — 
which we compute for n — 2 and n = 4 in Section Hm 
that depends on eccentricity, | | and |vfc|; this function 
can however not depend on the masses m„ m k of the 
particles, since the latter are decoupled from the veloci- 
ties due to the factorization of the space and momentum 
parts in the initial distribution ([7]). This means that the 
intrinsic dependence of the Fourier coefficient v n (pi) is 
on the velocity |v{|, rather than on the momentum pf. 
when comparing coefficients of different particle species, 
they take the same value at fixed velocity, i.e. as a func- 
tion of momentum they exhibit mass ordering. In order 
to emphasize the independence of the (qualitative) mass 
ordering of flow from the assumed initial spectrum, we 
shall present in Section ITVl results obtained with two dif- 
ferent choices for f(p). 

Once the anisotropics of the velocity distribution have 
been obtained, we can assume a specific form for the 
initial momentum distribution, and compute the flow co- 
efficients v n (pi). The denominator in definition (jlbl) is to 

our order of approximation given by replacing fa by f- ^ 
in the integral ©, which yields 

I d(pi SpT = Ni ^- (10) 

To ensure the internal consistency of our model, the 
average number of collisions per outgoing particle in the 
evolution should be small, say at most equal to 1. This 
number of collisions is quite naturally the integral over 
time of the collision rate, which for two beams of particle 



densities m, n k with velocities Vj, reads 

with the densities given by Eq. ([2]) and the relative veloc- 
ity by Eq. ([5])- One recognize here (minus) the loss term 
of the Boltzmann equation (QJ, integrated over p^. The 
calculation is thus very similar to that of the flow coef- 
ficients, with the only difference that there is no cosntpi 
term — which amounts to taking n = — and that there 
is a final integration over |pj| — which is the equivalent 
of computing the "integrated flow" . Fixing the average 
number of collisions to 1 will of course be equivalent to 
fixing the value of the cross section Od (or rather, of the 
dimensionless quantity a^/R), which allows us to present 
numerical values for i>2 (pi) and v^(jpi) in Section ITVl 

III. EXPLICIT CALCULATIONS 

In this Section we detail our computations of the 
anisotropy of the velocity distribution and of the mean 
number of collisions per particle, whose principle was pre- 
sented in the previous Section. The reader who is not in- 
terested in technical issues but is willing to trust us may 
skip directly to Section ITVl which only relies on Eqs. (j2"Tj) 
and [ggp . 

More specifically, we shall derive the anisotropy of the 
distribution in for particles of species i induced by col- 
lisions with particles of a different species k. We shall ar- 
bitrarily call the former "diffusing particles" and the lat- 
ter "scattering centers" . In the most general case where 
both types are massive, the relative velocity (O can be 
rewritten after some simple algebra as 

v lk = c^[l-p t p k cos(<p t -<p k )] 2 - (l-/?2)(l-/32), 

(12) 

where $ = |vj|/c, (3k = l v fc|/c Note that this expres- 
sion is simpler when one of the j3 is equal to 1, i.e. for 
massless particles. Accordingly the corresponding calcu- 
lations, which we present in Appendices[Bland[Cl become 
much more straightforward. 

As argued below Eq. ([9]), we need only consider the 
loss term of the collision integral, computed with the 
free-streaming solution ^ and integrated over x. The 
starting point of the calculation is thus the kernel 



C(t,pi,p k ) = - j d 2 xdO/f° ) ^,x,p fc )/^ 0) (i,x,p fc )'y lfc cr d 
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fi(Pi) fk(Pk) V [1 - costo-V*)] - (1 - ~ PI) 
( [l + e cos 2ipi] (3i - 2[ecos2{<pi + <p k ) + cos2((p l -<p k )]l3 t f3k + [l + e cos 2ip k ] 



, (13) 
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where we have performed the Gaussian integral over x and that over the "solid" angle G. 

As explained in Section [TT1 this kernel allows us to compute the number of collisions — by integrating it over 
time, ipi, ipki \Pk\ and |pj| — as well as the anisotropic flow coefficients v n (pi) — by multiplying with cosnifi before 
integrating over tpi, and leaving out the final integration over \p%\. 

The first integral, over time ranging from to oo, is trivial 



dtC(t,pi,pk) 



NtNkCTdVT- 
8tt 3 / 2 R 



■ fiipi) fk(pk) 

[l - ft&cosfo-^)] 2 " (1 - Pf)(l - 131) 



,1/2 



[l +ecos2y?i]/3? - 2 [e cos 2 (ipi + ip k ) + cos2((p l -(p k )]p i p k + [l + ecos2<p fc ]/3| 



(14) 



which then has to be multiplied by cos mpi, possibly with n = for the average number of collisions. 
A convenient approach to computing the angular integrals is to perform the change of variables 



(<Pi,(pk) -> (f 



o ) ° = « 



(15) 



which involves a Jacobian that can be reabsorbed in the limits of integration. The angular part to be integrated then 
reads 



cos n 



(0P + 5)y/ (1 - PA cos 2,5) 2 ~ (1 - P 2 ){1 - P\) 
yfas cos 20p + bs sin 2(f)P + eg 



(16) 



with (5-dependent coefficients as = [{Pf + Pi) cos 26 — 2p i p k ]e, bs = {Pi — Pf)e sin 26, and c$ = Pf + Pi — 2p i p k cos 26. 
The simplest is to perform first the integration over cb p , which gives 



cosn{cb p + 6) dcj) p 



- n y/as cos 2<j)P + bs sin 2</>p + cs 



Vl + e 

~2 

3ne 2 



K 



2c 



1 



1 + ey ^P 2 -2p l p k cos25 + Pl 



Pj -2p i p k cos28 + Pl cos 4(5 
(PI-2PA cos2<5 + /3 fc 2 )3/2 
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2 -fid 5 |i 3; e 2 ) 

Pf - Apfp k cos 26 + ZPlPl cos 46 - 4ft/3| cos 6<5 + g| cos { 
X (Pl-2p i p k cos26 + Plfl 2 



for n = 0, 
for rt = 2, 

for n = 4, 
(17) 



where is the complete elliptic integral of the first kind 
and 2-Fi is the Gaussian hypergeometric function|f| 

The next integral, over 6 G [— n,ir], is tedious, yet fea- 
sible analytically. A first step is to factorize the terms 
that come with half-integer powers in a rational fraction 



3 The hypergeometric functions in the second and third line can 
also be rewritten in terms of combinations of powers of e and of 
complete elliptic integrals K and E of ^2e/(l + e). The formu- 
lations given in Eq. Ill7jl are more compact, and allow one to see 
more easily the small-e behavior. Reciprocally, one might rewrite 
the expression for n = using the identity 



--K\ 



2e 
l + e 



' 1 3 
.44 



of cos 2 6: 



\l-p t p k cos 2(5) — (1 — Pf)(^ — Pi) 
(Pf-2p i p k cos2(5 + /3 2 )"+i 



with 



2P i p k 



Xi 



X3 



/ (1 — X\ COS 2 (5)(1 — X2 cos 2 6) 

V (1 - x 3 cos 2 8) n + 1 ' 
1 + p t p k - ^J(l-P 2 L ){1-P 2 k ) 



(Pi+Pkf 



2pip k 



%2 = 



l + PA + y/il-Pm-Pl) 



(P l + Pk) 2 

^PiPk XiX 2 



(Pt + Pk) 2 p t Pk 



(18) 



(19) 



Next, one performs the change of variable u = cos 2 6, 
which yields extra terms in the denominator under the 



5 



square root: the resulting integrals are of the form 



ity: 



duV n (u) 



I (1 — xiu)(l — X2II) 
u(l - u)(l - x 3 u) n+1 ' 



(20) 



with V n a polynomial of degree n, coming from the 
numerators in Eq. (fl7|) . Note that the integral over 
S G [— 7r,7r] gives 4 times the integral over u € [0,1], 
yet there comes a factor 5 from the change of variables, 
which partly compensates for it. 

In the case n = 0, the polynomial is a constant, so that 
we are left with a hyperelliptic integral. The latter can 
be expressed in terms of the Lauricella hypergeometric 
function of three variables , which is discussed in 
Appendix [A] Restoring the prefactors from Eqs. ([14"]) 
and (1171) , and dropping the minus sign which is irrelevant 
here, one finds 



NiN k a d ~ 



fi(Pi)fk(Pk)VT^eK\ 



F. 



W 2 



-1 -1 1 
~2> 2 



2c 



, 1; xx,x-z, x 3 



(21) 



To obtain the total number of collisions between particles 
of species i and k, one has to assume some shape for the 
initial spectra fi(pi), fk(j>k) and then to integrate over 
Pi = mcfii I y 1 — (3f and p k ■ Given expression (|21[) , in 
which the velocities enter through the variables x±, x-r. 
x 3 , an analytical calculation seems highly improbableO 
Nevertheless, one can use the fact that the Lauricella 
function in the second line actually only takes values be- 
tween 2/n and 1, when (3i and f5 k vary between and 
1, together with the normalization of fi and f k , to set 
bounds on the number of collisions at a given eccentric- 



2NiN k a d 

7T 3 / 2 i? 



2e 



1 



NiNkad 



2c 



Dividing by N t then yields the average number of col- 
lisions per particle of type i. If one wants to ensure 
that there is at most one collision per particle, then the 
first inequality translates into a bound on <j&. As the 
eccentricity-dependent part decreases with e from ir/2 
to 0, one can fix this bound to ensure less than one 
collision per particle over the whole eccentricity range: 

^r x = v^R/Nk. 

For the cases n = 2 and n — 4, a convenient trick is to 
rewrite V n {u) as a polynomial in (1 — x 3 u), as e.g. 



w = 4 



iPt+Pkf 



{l-x 3 uf 

(f3 t + l3k) 2 



- P k (l - x 3 u) 

which allows us to express the result of the integrals as 
sums of Lauricella functions of three variables only. In- 
cluding the prefactors from Eqs. (fl"4"]) and (fT7|). one finds 

Jdtd(p l dip k C(t,p. l ,p k )cosnip i = 

Nih{ Pl ) N k f k (p k )M n lC n (e)F n {Pi,p k ), (22a) 

with 



8R ' 



64i? ' 



(22b) 

Me) = V / T T ^ 2 Fi^,|2; £ ^ e (22c) 
/C 4 (e) = v / T T ^ 2 Fi^,^;3;e 2 )e 2 , (22d) 



and 



1 m+Pkf (3) a -1 -1 -1 \_ 2 (3 )/i ^i ± 1 

/3 2 [ 2 V 2 ' 2 ' 2 ' 2 '^i'^,^ P fc ^ D ^2' 2 ' 2 '2' ' 

(ft-^ (3) /l -1-13 \ 



xi,x 2 ,a;3 



(22e) 



4 Unless the integral over p fc and/or pt can be computed earlier in investigated this possibility, 

the calculation, before some of the other integrals. We have not 
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(Pi + ft) 4 (3)/l zl zl zl 
2 D W 2 ' 2 ' 2 



,l;xi,X2,x 3 -2ft(/3. 



n2 (3)/ 1 - 1 "I - 1 \ 

+ ft) I 2' ~2~' ~2~' ~2~' 1; ^i'^,^ ) 



+ Pi(ZPl-2pf)F l 



zi zi I . 

15 ^2' 2 ' 2 ' 2' ' 1 '" 



x 2 ,x 3 \-2Pi{fi l -p k YF i 



zi Zi 3 !. 

D ^2' 2 ' 2 ' 2' ' x * 2 ' 



3'3 



{Pi^Mi pwfi zi zi ! i. r 

2 y 2 ' 2 ' 2 ' 2' ' *' 



(22f) 



We shall exploit these lengthy formulae in the following 
Section, and give alternative derivations of the limiting 
cases pi = ft = 1 (resp. ft — 1, ft — 0) in Appendix iBl 
(resp.[C|. Before that, let us make two remarks. 

First, it is actually possible to derive shorter expres- 
sions, if instead of rewriting the polynomial V n (u) as we 
did here, one rather looks for its zeros, factorizes the 
polynomial, and inserts the factorized V n (u) under the 
square root, which yields n extra terms in the numera- 
tor of the rational fraction. The resulting integral can 
then be expressed as a single Lauricella hypergeometric 
function of 5 (resp. 7) variables for n = 2 (resp. n = 4). 
This is indeed more compact, yet much more expensive 
to compute, which explains our choice. 

The second remark is that Eqs. (|22ep - (|22f|) may not 
always be the most tractable expressions in practice, in 
particular when we shall be investigating the behavior of 
the J- n in the limiting cases ft <C 1 or ft = 1. As a matter 
of fact, they involve triply infinite series of variables that 
do not scale as a simple power of the velocities. It is then 
simpler to work at the level of Eq. (|16|) multiplied with 
the overall prefactor from the first line of Eq. (|14[) and 
with the relative velocity — performing e.g. the Taylor 
expansion in ft or considering the limit ft = 1 — , and 
then to integrate over 5. 



IV. RESULTS 

Let us now exploit the formulae (|2"Tj) - (l2"2")) derived in 
the previous Section. 



A. Eccentricity and velocity dependence 

In the previous Section, we have derived analytical 
formulae for the anisotropy of the velocity distribution, 
which integrated over |p^| yields the numerator in the 
definition (jTbJ) of v n (j>i), while the denominator is sim- 
ply given by Eq. (|TU| : 

v n (Pi) =A/"„/C„(e) J dp k p k N k fk(pk) Fn(Ph ft)- ( 23 ) 

Before we perform this last integration, which requires 
some ansatz for the initial momentum distribution, we 
can analyze some properties of the functions K, n and J- n . 



First, if one fixes the differential cross section to en- 
force at most one collision in average per particle, i.e. 
according to the discussion below (|2"T1) if o& = n^/nR/N^ 
with k < 1, then the prefactors M n simplify, and the 
whole dependence of v n (Pi) on N{, Nk, R cancels out. 

Next, the dependence on eccentricity is entirely given 
by the functions /C„, Eqs. p2c |) — (|22d |) ■ That is, v 2 (Pi) 
scales like e for small eccentricities, with less than 5% 
deviation from the linear scaling up to e = 0.75. v^(j>i) 
scales like e 2 for small eccentricities, with less than 5% 
deviation from the quadratic behavior up to e = 0.45. 
For larger eccentricities, K2 and /C4 fall off to reach for 
e = 1. 

An important fact is that the flow harmonics (|23|) 
depend on the particle velocity ft, rather than on its 
momentum pi. This clearly follows from the fact that 
the Boltzmann equation involves the velocities, not the 
momenta. As a consequence, particles with different 
masses but with the same velocity will have the same 
anisotropic flow w„(ft): when considered as a function 
of momentum, this gives rise to mass ordering of the 
flow coefficients, as is observed in experimental measure- 
ments For "slow" particles in the non-relativistic 
regime, ft ~ pi/mc, so that in our model we find that 
the flow coefficients should coincide at the same Pi/rrii 
for different particle species. This is noticeably what was 
found for slow particles emerging from an ideal expand- 
ing fluid [l2j], that is for particles that underwent many 
rescatterings. 

Let us discuss the dependence of F n on the velocity 
ft. At small ft and non- vanishing ft0 one finds the 
behaviors 

MPi,Pk) - y, Fi(PuPk) - (24) 

These scalings will not be modified in the integration over 
|pfe|. Since a small ft means a momentum pi ~ mcft, 
one thus finds that in the region where particles i are 



5 For /?£ = 0, which implies x\ = X2 = X3 = and corresponds to 
the scattering of particles of type i on fixed centers of type fc, the 
velocity-dependent parts in Eqs. l|22e|l - l|22fjl yield 1 for ft ^ 0. 
The limit ft = /3j. = is thus singular, yielding different results 
according to the order of the limits. This is rather obvious from 
the physical point of view, since it corresponds to both particles 
staying at rest. 
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non-relativistic, 

«n(Pi)ocp?, (25) 

which is the expected behavior when the particle momen- 
tum distribution is smooth at p = [10]. 

We show T% and —T4 (to take into account the minus 
sign of A/4) as functions of ft, ft in Fig- HI At fixed ft, 
one finds as expected the scalings (|2"4")1 at small ft — the 
factor 1/128 for T4 makes it difficult to display the quar- 
tic behavior, yet it is present. These smooth behaviors 
persist up to ft = ft, i.e. as long as the "diffusing par- 
ticles" i do not catch up with the "scattering centers" k 
emitted in the same direction. 

When ft > ft, the diffusing particles travel fast 
enough to probe the whole geometry of the distribution 
of scattering centers. There is thus for ft = ft a tran- 
sition from a regime where anisotropic flow is built up 
"locally" , involving only the diffusing particles that were 
close enough to the edge of the interaction region to find 
quickly a scattering center, to a "global" regime where 
every diffusing particle has a chance to scatter on an ini- 
tially distant center. As a consequence, the functions T2 
and J-4 vary sharply for ft > ft. 

F4 is especially interesting, for it changes sign. When 
the scattering centers barely move, ft <C ft, then the 
diffusing particles quickly probe the initial Gaussian dis- 
tribution ([7]). Now, a property of the latter is that 

(r 2 cos 2 ip r ) (r 4 cos4< f o r ) 3e 2 

P) = p) = 2T^' 

where (r, cp r ) denote the polar coordinates of a point in 
position space: x — r cos (p r , y = r sin tp r . One conven- 
tionally adds a global minus sign, to obtain a positive 
elliptic flow for e > [see Eq. (j8b[) ]. In turn, this means 
that one would naturally expect that a Gaussian distri- 
bution gives rise to a negative V4. This is exactly what 
comes out for —J- 4, which becomes negative when ft is 
significantly larger than ft. 



The flow harmonics at a given ft follow from integrat- 
ing F n over the whole ft range, Eq. (|23|) . Depending 
on the shape of the spectrum fk(pk), more weight is put 
on different ft regions, resulting in different values of 
v n (pi). This dependence on the spectrum will be partic- 
ularly marked for Vi(pi) due to the non-monotonic shape 
of -J" 4 (ft,ft) at fixed ft. 

Let us finally discuss the behavior of the flow coef- 
ficients in limiting regimes accessible in our approach. 
Corresponding direct computations, which are signifi- 
cantly simpler than going through Section Mil and then 
taking the relevant limits, can be found in Appendices IB"1 
and ED 

The first limiting regime corresponds to the massless 
gas limit, in which both ft and ft are uniformly equal 
to 1 (Appendix [Bj . The flow coefficients show an univer- 
sal behavior, i.e. v n is independent of the assumed initial 
spectrum. When the cross section is fixed so that parti- 
cles rescatter at most once in average (and exactly once 
in average for central collisions), we find U2A — 0.083 
with less than 5% deviation up to e = 0.75. 

The other limiting regime corresponds to that of a 
Lorentz gas of massless diffusing particles scattering on 
infinitely massive fixed centers, that is ft = 1, ft = 
(Appendix [Cj . In that case the diffusing particles probe 
most cleanly the Gaussian distribution of scattering cen- 
ters, and accordingly one finds a positive V2 and a neg- 
ative 1)4, again universally independent of the initial 
momentum distribution, in perfect agreement with the 
initial geometry. If one requires a single rescattering 
per particle in average for e = 0, then v-i ~ 0.25e, 
V4 ~ — 0.094e 2 , that is Ui/wf — —1.5, across a wide e 
range. 



B. Momentum dependence 

In order to present the flow coefficients as a function 
of momentum, as is customary, we now have to assume 
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some initial spectrum fk{pk) to insert in the integral 
We choose two different distributions: 

• A thermal-like spectrum 



hipk) = * e-y/X+<*' T , (26) 

1 (mfcc + 1 ) 

where we fix the inverse slope parameter T so as 
to give pions an average (transverse) momentum 
(p n ) = 420 MeV/c. 

• A QCD-inspired power-law spectrum, regulated at 
the origin 



fk(Pk) - — 2 

kpo i + (p±y 



(27) 



where the characteristic momentum scale po is 
fixed such that the average momentum per particle 
equals 420 MeV/c. 

Integrating Eq. (|2U)) with these spectra necessitates a nu- 
merical approach. 

To emphasize the mass ordering which we discussed 
above, we consider a mixture of three components with 
different masses (and different yields, to mimic a realis- 
tic situation, although we did not attempt to take ex- 
perimental values): 80% pions, 12.5% kaons and 7.5% 
protons. Accordingly, the mixture is described by a sys- 
tem of three coupled Boltzmann equations Q , which all 
in all involve 6 different collision terms (ir-n, 7T-K, 7r-p, 
K-K, K-p and p-p) that have to be taken into account 
carefully to ensure a fixed mean number — which will 
be taken equal to 1 for e = — of collisions per parti- 
cle. In turn, the flow harmonic v„(pi) for a given species 
involves three different terms. 

To keep the number of free parameters to a minimum, 
we assume the same differential cross section ad for all 
these collisions, and the same inverse slope parameter T 



or momentum scale po for all species. Note that this 
implies that the average momentum of heavier parti- 
cles is not the same for both choices of spectra: for the 
power-law spectrum, all species have a mean momentum 
of 420 MeV/c, while this only holds for pions for the 
thermal-like spectrum. 

We display our results for v 2 {pi) in Fig. [5] for both 
choices of initial spectrum, in the case e = 0.1. As antici- 
pated, there is for both initial momentum distributions a 
clear mass-ordering of v 2 {jpi), which is thus no hallmark 
of the presence of a thermal system. 

Another striking feature is the size of the elliptic flow, 
of the order of a few percent, i.e. v 2 /e of the order of 0.1- 
0.2 for both spectra. Note that the magnitude of v 2 for a 
given species depends on its abundance, since we enforce 
the condition of an average single collision per particle on 
all particles, not for each species individually. Then, the 
more abundant a species, the more collisions it undergoes 
and thus the larger its elliptic flow. For that reason, the 
absolute value of v 2 for a given initial momentum distri- 
bution is not very meaningful in our eyes — it depends 
on the particle abundance, and additionally it is directly 
proportional to <7d, thus determined by our requiring a 
given average number of rescatterings per particle, which 
is an arbitrary choice of ours. 

Eventually, one can recognize the quadratic growth at 
small pi — or rather, at small Pi/rrii — for kaons and 
protons. For pions, the non-relativistic regime where this 
scaling holds is not visible. 

As explained in the discussion of J-4, we expect a richer 
variety of behaviors for the fourth harmonic v^{jpi). Ac- 
cordingly, we present more plots. In Fig.[3]we show Uj(pi) 
for pions, kaons and protons for both initial momentum 
distributions: thermal- like on the left panel and power- 
law on the right panel. Figure [4] then shows vn(j>i) for a 
mixture of 88% pions, 8% A and 4% D-mesons, with a 
thermal-like initial spectrum, so as to have even heavier 
particles, with smaller at a fixed pk- 
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p T (GeV/c) Pj (GeV/c) 

FIG. 3. Hexadecupolar flow v&{pi) of pions (full line), kaons (dot-dashed line) and protons (dashed line) assuming a thermal-like 
initial spectrum (|26p (left) or a power-law spectrum (|27|l (right). 
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FIG. 4. Hexadecupolar flow V4,(pi) of pions (full line), A (dot- 
dashed line) and D-mesons (dashed line), assuming a thermal- 
like initial spectrum (1261 1, 



sponds to where JF 2 is largest, which agrees with the 
larger v^ipi) found with the power-law spectrum. 

When one replaces kaons and protons by smaller 
amounts of the significantly heavier A and D, there are 
more "slow" scattering centers in the system, i.e. one 
again puts more emphasis on the negative —IF 4 regions. 
As a result, the pions V4(p 7r ) for a thermal- like spec- 
trum now takes negative in a small region, before 
turning positive again (Fig. 2]). A small region in p v 
amounts to much larger regions in p\ or pd, so that V4 
for those species remains negative over most of the dis- 
played range Q 

For the sake of completeness, let us mention the ratio 
V4(pi)/v2(pi) 2 ■ We do not show our results for this ratio 
because we do not give them too much worth: in our 
model V4(pi) / V2(pi) 2 obviously scales as 1/ay and might 
not be reliable in the limit of a low number of collisions 
we are considering, as we have not estimated the size of 
the corrections to our results. 



Again, v^ipi) exhibits mass ordering in the low momen- 
tum, non-relativistic region, irrespective of the choice of 
spectrum. A further consequence of the dependence of 
«4 on ^ rather than on pi is mostly visible on the right 
panel of Fig. [3l although present on all our plots of ellip- 
tic and hexadecupolar flow. Here, one clearly sees that 
the proton and kaon V4 have exactly the same shape as 
the pion one, although dilated along the momentum axis; 
in particular, they all reach the same minimum value. 

The most striking feature when comparing both pan- 
els of Fig. [3] is obviously the change of sign of V4,(pi) ac- 
cording to the choice of spectrum. As explained above, 
this comes from weighting the various velocity regions of 
J"4(/3i, (3k) differently in the integral over |pfc|: the power- 
law spectrum weights the low-velocity region more than 
the thermal-like distribution, so that it emphasizes re- 
gions of negative — J-4 more, which explains the negative 
Vi(pi). Note that these regions of negative —J- 4 corre- 



V. DISCUSSION 

Before we start discussing our results, let us summa- 
rize the salient features of our study. We have solved the 
Boltzmann equation in the limit of a (very) small aver- 
age number of rescatterings per outgoing particle, which 
allows us to compute the anisotropic flow harmonics in 
that limit. We then find that sizable elliptic flow is gen- 
erated already after a single rescattering and that V2(pi) 
exhibits mass ordering. Additionally, we find that the 
fourth harmonic tti(pj) is quite sensitive to the initial 
momentum distribution of particles. 



6 Interestingly, the STAR Collaboration has reported values of V4 
for antiprotons and A which are compatible with or even neg- 
ative in the low momentum region [ill . Fig. 7] . 
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Our approach can only be considered as a toy model 
for the phenomenology of nucleus-nucleus collisions, so 
that one should only retain the qualitative aspects, and 
not attach too much worth to the quantitative values 
we obtain. A first reason is that we have deliberately 
reduced the numbers of ingredients and parameters to 
a minimum: Gaussian initial shape, hard-spheres-type 
universal isotropic cross section, and initial spectra de- 
pending on one parameter only. The second reason is 
that we have restricted ourselves to the leading term in 
an expansion in powers of the cross section of the solu- 
tion of the Boltzmann equation. If one wants to evaluate 
the quantitative uncertainties on our results, one has to 
explicitly compute the next term in the expansion, or at 
least to estimate it. This can be done, at least numeri- 
cally, yet we have not attempted to do it. 

Nonetheless, we feel confident that the qualitative 
trends we have found are quite robust, and must also be 
present in more precise studies. In particular, the mass 
ordering of v n (j>T) follows from the fact that the Boltz- 
mann equation framework only involves the velocities of 
particles, and thus its origin appears to be intrinsic to 
the kinetic equation. 

Let us eventually discuss possible generalizations of the 
present work. First, one may try different shapes for the 
initial distribution in position space, in particular more 
realistic, finite shapes. 

A desirable improvement that might become feasible 
would be the possibility to perform both angular inte- 
grals first, so as to isolate the time evolution of the flow 
harmonics — as is the case for the Lorentz gas, see Ap- 
pendix [UJ This would help implement our study as a 
model for two rather different cases: first, as a descrip- 
tion of the very first steps of the expansion, which would 
provide initial conditions for a subsequent hydrodynamic 
calculation. Second, one may possibly describe momen- 
tum anisotropy at high momentum: since we only con- 
sider the loss term of the (elastic) collision term, our ap- 
proach is not far from constituting a phenomenological 
description of the loss of partons with a given high mo- 
mentum due to their branching in two lower momentum 
partons. 

Another possible generalization consists in introducing 
some dependences in the cross section: either on the col- 
lision type, or on the center-of-mass energy of the rescat- 
terings, or allowing for anisotropic cross sections. 

Finally, taking care of the next order (s) in the expan- 
sion will be needed, if only to estimate the uncertainty 
on the results. 
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Appendix A: Lauricella's multivariate 
hypergeometric function Fd 

In this Appendix we present some basic properties of 

(3) 

the function with the help of which we express our 
analytical results for the anisotropy of the velocity dis- 
tribution. 

The Lauricella function of the fourth kind of m € N + 
variables Z\,...,z m is defined by (l3j 



7 (m) 

D 



E 

fei,...,fc m =0 



(a, b OT ,c; zi, 

i a )k 1 +---+k m (h)k 1 



k 1 



(c)fc 1+ ... + fc m &i! 

where the Pochhammer symbol (a)k is defined by 

T(a + k) 



(Al) 



{a) k 



T{a) 



< 1. 



The series converges absolutely for \zi\ < 1, . . . , \z„ 
For m — 1 , it reduces the Gauss hypergeometric function 
2-Fi(a, b, c, z), and for m — 2, to the Appell function F\. 

If the real parts of a and c—a are positive, the function 
admits the integral representation 13 1 

^1) • • • ) jCJ Z± , . . • , Z m ) — 

r(c) f 1 u-^—'du 
r(o) r( c - a) J (i - Zl u) b ^ ■■■(!- z m u) b ™ ■ K 1 

This representation allows the analytical continuation of 
Fp 1 ^ to the cut complex z^-planes deprived of the half- 
line from 1 to infinity along the positive real axis. 
If Re (c — a — bi — • ■ ■ — b m ) > one finds 



F { j^\a, & m ,c;l,...,l) = 

r(c) r(c — a — b\ 



b m ) 



r(c — a) r(c — b\ 



b m ) 



(A3) 



Appendix B: Massless diffusing particles and 
scattering centers 

In this Appendix we calculate the anisotropics of the 
velocity distribution and the average number of collisions 
in the case where both colliding particle species are mass- 
less, which leads to = 0k = 1. We assume that the 
particles remain distinguishable. 

In that case, the relative velocity (fT2"j) reduces to 



Vlk 



c[l — COs((fi 



2c sin 



2<Pi~ Vk 



(Bl) 



and the kernel f| 13 j) becomes, after performing the change 
of angular variables (TT5)) , 



C(t,Pi,p k ) 



NiN k a d cVl - e 2 - - 

- JiKPi) JkiPk) 



47T 2 i? 2 



x sin (5 exp 



cH 2 
R 2 



(l - ecos2</. p ) sin 2 6 



(B2) 
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which simplifies significantly the ensuing calculations and 
the final formulae compared to the case of arbitrary ve- 
locities. 

As above, we need to integrate the kernel (|B2[) multi- 
plied by cos n(4> p + 5) over both angles </> p and S and over 
time to obtain the average number of collisions (n = 0) 
and the flow coefficients. 

The integral over time is trivial and yields 



N z N k a d Vl - < 



8it 3 / 2 R 



2 ?, N r, N | sinS\ cosn((f) p + 5) 

JiijPi) Jk(Pk) 7= = , = — • 

yl — ecos20P 



The angular part of this expression should be compared 
with Eq. (|T6|. 

Integrating over S € [—it, it] is also straightforward and 
gives 



nl 1 i\ -ill r, Ji{Pi) Jk\Pk) 

2{n z — \)-n A l z R 



VI - ecos20P 



(B3) 



There remains the integral over <j) p € [— it, it]. One 
finds for n = 



2N t N k a d 2 



ir 3 / 2 R 



fiiPi) fkiPh) Vl - e # 



2e 



(B4) 



where we have dropped the minus sign, which comes from 
our integrating the loss term of the Boltzmann equation, 
and is irrelevant for the computation of the mean number 
of collisions. Now, if one takes ft = ft^ = 1 in Eq. ([2~T]) . 
then x\ — xi — x^ — 1 and the value of Lauricella func- 
tion in the second line is given by Eq. (|A3|) with a = A, 



bi 



6, = i 



1, i.e. it equals 2/tt, so that 



Eq. (|B4|) does represent the limit of Eq. ([21]) . 

Integrating Eq. (|B4I) over |pj| and |pfe| is straightfor- 
ward and one finds that the total number of collisions 
between particles i and k is 



N r , 



2N t N k a d 

7T 3 / 2 i? 



s/T~tK\ 



2e 
1 + e 



If one requires that each "diffusing particle" i undergo 
in average one collision, so that N co \\ = N%, for the most 
central collisions, where the eccentricity-dependent term 
equals it/ 2, then this fixes the differential cross section 
to 



-R. 



(B5) 



Let us now turn to the anisotropic flow harmonics. The 
integral of Eq. (|B3[) over 4> p for n — 2 or n = 4 reads after 
some algebra 

|^/,Wih)^.F,(y;2; e j e , (B6a) 
^/ i(K )/ lW VI^ 2Fl (^ !3ie ^. ( B6b) 



One recognizes the eccentricity-dependent (|22c[) (|22dl) . 
Furthermore using Eq. (|A3|) with a — h, b-y =62 = — h, 
c = 1 and &3 = j, — A or — |, one finds ^(l, 1) = 2/37T 



and 74(1,1) = 2/157T, so that Eqs. (|22J reduce to 
Eqs. (p| for ft = k = 1. 

Again, the latter can be integrated over |pfc|, divided 
by the denominator (JTU]) of the flow coefficient definition, 
which yields the flow harmonics 



V 2 (Pi) 



N k a d 



12v^i? 
N k a d 

leov^i? 



V / T^ 2 ^fil;3;e 2 ) e 



4' 4" 



Note that these values are actually independent of the 
momentum, and in particular do not vanish when pi goes 
to 0! This comes from the facts that in our model the 
harmonics v n depend primarily on velocity, rather than 
on momentum, and that massless particles always travel 
with the velocity of light, irrespective of their momentum. 

If the cross section is fixed to its value (|F35|) to ensure 
one collision per particle in the most central collisions, 
these flow harmonics becomes 



l>2 



Vi{Pi 



12 \4 4 

'5 7 
I'! 5 



— y/l-e 2 2 F 1 ( 
160 V \ 



3;6 2 W 2 . 



This gives a rather large U2/e — 0.08 for e up to 0.75, 
and a ratio v^/v\ that equals 0.9 for e = and slowly 
increases quadratically with e. 



Appendix C: Massless diffusing particles scattering 
on fixed centers 



We now consider the case of a Lorentz gas of massless 
diffusing particles (ft = 1) scattering on fixed, infinitely 
massive centers (ft. = 0). 

In that case, the relative velocity (| 12[) is Vi k — c, in- 
dependent of the relative angle tpi — tp k , which is quite 
natural as tp k is undefined. In turn, the kernel (|13[) sim- 
plifies to 



C(t,pi,p k ) = 



8n 2 R 2 
c 2 t 2 
AR 2 



■ fiiPi) fk(Pk) 



x exp 



1 + e cos 2(f i 



independent of ip k , over which one can at once integrate, 
which amounts to multiplying the kernel with 27r. The 
resulting expression, multiplied with cos rupi , is easily in- 
tegrated over ifi, yielding 



(-1) 



l/2- 



_! NiN k a d cVl - e 2 7 , N 7 , N 
— fiiPi) Jk(Pk) 



x e 



2i? 2 

-c 2 t 2 /4R 2 



cH 2 
\4R 2 ' 



(CI) 
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where we have assumed that n is even, and where I n /2 
denotes the modified Bcssel function of the first kind of 
order n/2. Note that the factor (— 1)"/ 2 + 1 comes in part 
from the global minus sign due to our considering the 
loss term of the Boltzmann equation, and in part from 
our writing the Bessel function with a positive argument. 

Now, integrating Eq. (jClj) over time and over |p^| and 
dividing by the denominator (|10[) yields the flow har- 
monic v n (pi), that is the large-time limit of v n (pi,t). The 
integral over time of the second line of Eq. (|C1|) gives [H], 
formula 2.15.3.2] 

c 2«+i(f)! £ 2F \— '2 +1,C J' (C2) 

with (2n - 1)!! = 1 • 3 . . . (2n - 1) if n > 1, 1 if n = 0. 

Taking into account the prefactor from Eq. (|C1[) and 
performing the trivial integrals over |p/-| and |pj| gives 
with n = the total number of collisions 



NiN k <T d yfif 

2R 

N l N k a d 



Vl-^aFil T , T ;l;e' 



1 3 

4' 4 



2e 



which is indeed the limit of Eq. (f2"Tj) for j3 k — — i.e. 
x\ = X2 — £3 — — after integrating out the initial 
spectra. To ensure one collision per diffusing particle for 
e = 0, and thereby less than one collision over the whole 
eccentricity range, the cross section should be 



1 coll 



R. 



(C3) 



Using Eq. (|C2j) with n = 2 or n = 4 and restoring the 
necessary prefactor from Eq. (jCip . one derives the flow 
harmonics 



N k a d y/TT 



Vi{pi) = - 



8R 
3N k a d ^/jr 
64i? 



4' 4' 
5 7 



(CM) 



Vl-e^Fil j, i; 3;e^)e^, (C5) 



which is what Eqs. (|22|) give for f3 k — 0, ^, 7^ 0. Replac- 
ing ad by its value (|C3|) gives 



'-'2 



5 7 



That is, the diffusion of massless particles on fixed scat- 
tering centers leads to a large, momentum independent 
v 2 /e = 0.25 for e up to 0.75, while is negative, which 
directly reflects the geometry of the distribution, as dis- 
cussed in Section TlVBI This leads to a ratio v±lv\ = 
-1.5. 

If we do not integrate Eq. (|C1I) over time, but perform 
the integration over |pfe| and the division by Nifipi), we 
obtain the time derivative dv n /dt: 



dv n = , l y l /2 + l N kPdC 



dt 



s/T 



R? 



J <4^ £ J' 



This derivative is interesting, for it provides us at once 
with the small-time evolution of the flow coefficients, 
valid for t < 2R/c: 



dv n , 

dt ~ [ ~ 



m\R 2 \ 4R J 



i.e. after integration 

v n (t) oc (-lp 2 +H n+1 . 

One thus recovers the behavior found in Ref. [j| and con- 
firmed in Monte Carlo simulations [15]. We were un- 
fortunately not able to derive this result analytically in 
the more general of collisions of massive particles, by 
performing both angular integrals before that over time 
without making some additional assumption (as e.g. an 
expansion in small eccentricities). 
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